function [x,output]=agentsolverho_4agents(xguess,index,ewsb,waut,bind,nobind);
global P beta n m S state agent gammagrid  Y scalefactor lowerbound rho table
objval = [];
output=[];
options = optimset('LargeScale','off','Display','iter');
[x,fval,exitflag]=fmincon(@objfun,xguess,[],[],[],[],[],[],@constrfun,options);
output(6)=exitflag;
function f = objfun(x);
i=index(1);
j=index(2);
consu(1:agent-1) = x(1:agent-1);
consu(agent)=Y(j)-sum(x(1:agent-1));
wint=interpolation_4agents(ewsb, index,x./Y(j));
for g=1:agent
output(g)=utility(rho,consu(g))+beta*wint(g);
end
objval = gammagrid(i,1)*output(1) + gammagrid(i,2)*output(2) + gammagrid(i,3)*output(3) + gammagrid(i,4)*output(4);
f=-objval;
end
function [cineq,ceq]=constrfun(x);
i=index(1);
j=index(2);

consu(1:agent-1) = x(1:agent-1);
consu(agent)=Y(j)-sum(x(1:agent-1));

ceq=[];
for k=1:length(nobind)
ceq = [ceq; gammagrid(i,nobind(end))/gammagrid(i,nobind(k))-consu(nobind(end))/consu(nobind(k))];
end
for k=1:length(bind)
ceq = [ceq; ((utility(rho,S(j,bind(k)))+beta*P(j,:)*waut(:,bind(k))) - output(bind(k)))];
end
cineq=[];
end
end